#include "field.h"

void  field_set_pressure(Field           *f_P,
						 state_equation   eq,
						 const Field     *f_rho,
						 const double     T)
{
    CHECK_FIELD(f_P);
    CHECK_FIELD(f_rho);
    CHECK_FIELDS(f_P,f_rho);
    assert(NULL!=eq);
    
    {
        double **       P   = f_P->array;
        double ** const rho = f_rho->array;
        const int       nx  = f_P->nx;
        const int       ny  = f_rho->ny;
        int i,j;
        
        for(j=0;j<ny;++j)
        {
            for(i=0;i<nx;++i)
            {
                P[j][i] = eq( rho[j][i], T );
            }
        }
    }
    
}

void field_set_viscosity(Field *f_eta, 
						 state_equation visco,
						 const Field *f_rho,
						 const double T)
{
	CHECK_FIELD(f_eta);
	CHECK_FIELD(f_rho);
	CHECK_FIELDS(f_eta, f_rho);
	{
		double ** eta1 = ARRAY_OF(f_eta);
		const double **rho = CONST_ARRAY_OF(f_rho);
		const int nx = f_eta->nx, ny = f_eta->ny; 
		int i,j; 
		
		for(j=0;j<ny;++j)
        {
            for(i=0;i<nx;++i)
            {
                eta1[j][i] = visco( rho[j][i], T );
            }
        }
	}
}

void field_set_second_viscosity(Field *f_zeta,
						   state_equation second_visco,
						   const Field *f_rho,
						   const double T)
{
	CHECK_FIELD(f_zeta);
	CHECK_FIELD(f_rho);
	CHECK_FIELDS(f_zeta, f_rho);
	{
		double ** zeta1 = ARRAY_OF(f_zeta);
		const double **rho = CONST_ARRAY_OF(f_rho);
		const int nx = f_zeta->nx, ny = f_zeta->ny; 
		int i,j; 
		
		for(j=0;j<ny;++j)
        {
            for(i=0;i<nx;++i)
            {
                zeta1[j][i] = second_visco( rho[j][i], T );
            }
        }
	}
}

void field_compute_momenta(Field       *f_Ju,
                           Field       *f_Jv,
                           const Field *f_rho,
                           const Field *f_u,
                           const Field *f_v)
{
    CHECK_FIELD(f_Ju);
    CHECK_FIELD(f_Jv);
    CHECK_FIELD(f_rho);
    CHECK_FIELD(f_u);
    CHECK_FIELD(f_v);

    CHECK_FIELDS(f_Ju,f_Jv);
    CHECK_FIELDS(f_Ju,f_rho);
    CHECK_FIELDS(f_rho,f_u);
    CHECK_FIELDS(f_rho,f_v);

    {
        double **       Ju  = ARRAY_OF(f_Ju);
        double **       Jv  = ARRAY_OF(f_Jv);
        const double ** rho = CONST_ARRAY_OF(f_rho);
        const double ** u   = CONST_ARRAY_OF(f_u);
        const double ** v   = CONST_ARRAY_OF(f_v);
        const int       nx  = f_rho->nx;
        const int       ny  = f_rho->ny;
		const int       nx1 = nx-1; 
		const int       ny1 = ny-1; 
        int i,j;
        for(j=0;j<ny;++j)
        {
           // const int jm = (j <= 0)   ? ny1 : j-1;
            const int jp = (j >= ny1) ? 0   : j+1;
            
            for(i=0;i<nx;++i)
            {
                //const int im = (i <= 0)   ? nx1 : i-1;
                const int ip = (i >= nx1) ? 0   : i+1;
		        const double rho0 = 0.5*(rho[j][i]+rho[jp][i]);
				const double rho1 = 0.5*(rho[j][ip]+rho[j][i]); 
                Ju[j][i] =  rho1 * u[j][i];
                Jv[j][i] =  rho0 * v[j][i];
            }
        }
        
    }
}

void field_compute_velocities(Field       *f_u,
                              Field       *f_v,
                              const Field *f_rho,
                              const Field *f_Ju,
                              const Field *f_Jv)
{
    CHECK_FIELD(f_Ju);
    CHECK_FIELD(f_Jv);
    CHECK_FIELD(f_rho);
    CHECK_FIELD(f_u);
    CHECK_FIELD(f_v);
    
    CHECK_FIELDS(f_Ju,f_Jv);
    CHECK_FIELDS(f_Ju,f_rho);
    CHECK_FIELDS(f_rho,f_u);
    CHECK_FIELDS(f_rho,f_v);
    
    {
        const double ** Ju  = CONST_ARRAY_OF(f_Ju);
        const double ** Jv  = CONST_ARRAY_OF(f_Jv);
        const double ** rho = CONST_ARRAY_OF(f_rho);
        double **       u   = ARRAY_OF(f_u);
        double **       v   = ARRAY_OF(f_v);
        const int       nx  = f_rho->nx;
        const int       ny  = f_rho->ny;
		const int       nx1 = nx-1; 
		const int       ny1 = ny-1; 
        int i,j;
        for(j=0;j<ny;++j)
        {
           const int jm = (j <= 0)   ? ny1 : j-1;
            const int jp = (j >= ny1) ? 0   : j+1;
            
            for(i=0;i<nx;++i)
            {
                const int im = (i <= 0)   ? nx1 : i-1;
                const int ip = (i >= nx1) ? 0   : i+1;
				
		        const double rho0 = 0.5*(rho[j][i]+rho[jp][i]);
				const double rho1 = 0.5*(rho[j][ip]+rho[j][i]); 
                u[j][i] = 0.5*(Ju[jp][i]+Ju[jm][i])/rho1;
                v[j][i] = 0.5*(Jv[j][ip]+Jv[j][im])/rho0;
            }
        }
        
    }

}

void field_compute_increment(Field *f_delta_Ju,
							 const Field *f_dvg,
							 const Field *f_grad,
							 const Field *f_rx)
{
	CHECK_FIELD(f_delta_Ju);
	CHECK_FIELD(f_dvg);
	CHECK_FIELD(f_grad);
	CHECK_FIELD(f_rx);
	CHECK_FIELDS(f_delta_Ju, f_dvg);
	CHECK_FIELDS(f_delta_Ju, f_grad);
	CHECK_FIELDS(f_delta_Ju, f_rx);
	{
		double       **delta_Ju  = ARRAY_OF(f_delta_Ju);
		const double **dvg       = CONST_ARRAY_OF(f_dvg);
		const double **grad_P    = CONST_ARRAY_OF(f_grad);
		const double **rx        = CONST_ARRAY_OF(f_rx); 
		
		const int nx = f_delta_Ju->nx; 
		const int ny = f_delta_Ju->ny; 
		int i,j; 
		
		for( j=0; j < ny; ++j )
		{
			for( i=0; i < nx; ++i )
			{
				delta_Ju[j][i] = dvg[j][i]-grad_P[j][i]-rx[j][i];
			}
		}
	}
}

void field_compute_increment_1(Field *f_delta_Ju1,
							   const Field *f_dvg,
							   const Field *f_grad,
							   const Field *f_laplacien,
							   const Field *f_rx,
							   const double eta)
{
	CHECK_FIELD(f_delta_Ju1);
	CHECK_FIELD(f_dvg);
	CHECK_FIELD(f_grad);
	CHECK_FIELD(f_rx);
	CHECK_FIELD(f_laplacien);
	CHECK_FIELDS(f_delta_Ju1, f_dvg);
	CHECK_FIELDS(f_delta_Ju1, f_grad);
	CHECK_FIELDS(f_delta_Ju1, f_rx);
	CHECK_FIELDS(f_delta_Ju1, f_laplacien);
	{
		double       **delta_Ju1 = ARRAY_OF(f_delta_Ju1);
		const double **dvg       = CONST_ARRAY_OF(f_dvg);
		const double **grad      = CONST_ARRAY_OF(f_grad);
		const double **Lrhu      = CONST_ARRAY_OF(f_laplacien);
		const double **rx        = CONST_ARRAY_OF(f_rx);
		
		const int nx = f_delta_Ju1->nx, ny = f_delta_Ju1->ny; 
		int i,j; 
		
		for( j=0; j < ny; ++j )
		{
			for( i=0; i < nx; ++i )
			{
				delta_Ju1[j][i] = eta*Lrhu[j][i]+dvg[j][i]-grad[j][i]-rx[j][i];
			}
		}
		
	}
	
}